2-1. Dependency Measures
library(reshape2) # melt function
library(ggplot2) # ggplot function
library(pcaPP) # Fast Kendall function
library(energy) # Distance Correlation
library(Hmisc) # Hoeffding's D measure
## Loading required package: lattice
## Loading required package: survival
## Loading required package: Formula
##
## Attaching package: 'Hmisc'
## The following objects are masked from 'package:dplyr':
##
## src, summarize
## The following object is masked from 'package:SeuratObject':
##
## Key
## The following object is masked from 'package:Seurat':
##
## Key
## The following objects are masked from 'package:base':
##
## format.pval, units
library(zebu) # Normalized Mutual Information
# library(minerva) # Maximum Information Coefficient
library(XICOR) # Chatterjee's Coefficient
# library(dHSIC) # Hilbert Schmidt Independence Criterion
library(VineCopula) # Blomqvist's Beta
make_cormat <- function(dat_mat){
matrix_dat <- matrix(nrow = ncol(dat_mat), ncol = ncol(dat_mat))
cor_mat_list <- list()
basic_cor <- c("pearson", "spearman")
# find each of the correlation matrices with Pearson, Spearman, Kendall Correlation Coefficients
for (i in 1:2){
print(i)
cor_mat <- stats::cor(dat_mat, method = basic_cor[i])
cor_mat[upper.tri(cor_mat, diag = T)] <- NA
cor_mat_list[[i]] <- cor_mat
}
# functions that take matrix or data.frame as input
no_loop_function <- c(pcaPP::cor.fk, Hmisc::hoeffd,
VineCopula::BetaMatrix)
for (i in 3:5){
print(i)
fun <- no_loop_function[[i-2]]
cor_mat <- fun(dat_mat)
if (i == 4){ # Hoeffding's D
cor_mat <- cor_mat$D
}
cor_mat[upper.tri(cor_mat, diag = T)] <- NA
cor_mat_list[[i]] <- cor_mat
}
# functions that take two variables as input to calculate correlations.
need_loop <- c(zebu::lassie, energy::dcor2d, XICOR::calculateXI)
for (i in 6:8){
print(i)
fun <- need_loop[[i-5]]
cor_mat <- matrix(nrow = ncol(dat_mat),
ncol = ncol(dat_mat))
for (j in 2:ncol(dat_mat)){
for (k in 1:(j-1)){
if (i == 6){
cor_mat[j, k] <- fun(cbind(dat_mat[, j], dat_mat[, k]), continuous=c(1,2), breaks = 6, measure = "npmi")$global
} else {
cor_mat[j, k] <- fun(as.numeric(dat_mat[, j]),
as.numeric(dat_mat[, k]))
}
}
}
cor_mat[upper.tri(cor_mat, diag = T)] <- NA
cor_mat_list[[i]] <- cor_mat
}
return(cor_mat_list)
}
draw_heatmap <- function(cor_mat){
len <- 6
melted_cormat <- melt(cor_mat)
melted_cormat <- melted_cormat[!is.na(melted_cormat$value),]
break_vec <- round(as.numeric(quantile(melted_cormat$value,
probs = seq(0, 1, length.out = len),
na.rm = T)),
4)
break_vec[1] <- break_vec[1]-1
break_vec[len] <- break_vec[len]+1
melted_cormat$value <- cut(melted_cormat$value, breaks = break_vec)
heatmap_color <- unique(melted_cormat$value)
heatmap <- ggplot(data = melted_cormat, aes(x = Var2, y = Var1, fill = value))+
geom_tile(colour = "Black") +
ggplot2::scale_fill_manual(breaks = sort(heatmap_color),
values = rev(scales::viridis_pal(begin = 0, end = 1)
(length(heatmap_color)))) +
theme_bw() + # make the background white
theme(panel.border = element_blank(), panel.grid.major = element_blank(),
panel.grid.minor = element_blank(), axis.ticks = element_blank(),
# erase tick marks and labels
axis.text.x = element_blank(), axis.text.y = element_blank())
return (heatmap)
}
make_cor_heatmap <- function(dat_mat, cor_mat_list){
fun_lable <- c("Pearson's Correlation", "Spearman's Correlation", "Kendall's Correlation",
"Hoeffding's D", "Blomqvist's Beta", "NMI",
"Distance Correlation", "XI Correlation")
cor_heatmap_list <- list()
cor_abs_heatmap_list <- list()
# make correlation matrices
#cor_mat_list <- make_cormat(dat_mat)
for (i in 1:8){
print(i)
if (i == 5){
next
}
cor_mat <- abs(cor_mat_list[[i]])
# get heatmaps
cor_heatmap <- draw_heatmap(cor_mat)
# ggplot labels
ggplot_labs <- labs(title = paste("Heatmap of", fun_lable[i]),
x = "",
y = "",
fill = "Coefficient") # change the title and legend label
cor_heatmap_list[[i]] <- cor_heatmap + ggplot_labs
if (i %in% c(1,2,3,4,6)){
cor_abs_mat <- abs(cor_mat_list[[i]])
cor_abs_heatmap <- draw_heatmap(cor_abs_mat)
ggplot_abs_labs <- labs(title = paste("Abs Heatmap of", fun_lable[i]),
x = "", # change the title and legend label
y = "",
fill = "Coefficient")
cor_abs_heatmap_list[[i]] <- cor_abs_heatmap + ggplot_abs_labs
} else {
ggplot_abs_labs <- labs(title = paste("Abs Heatmap of", fun_lable[i]),
subtitle = "Equivalent to Non-Abs Heatmap",
x = "", # change the title and legend label
y = "",
fill = "Coefficient")
cor_abs_heatmap_list[[i]] <- cor_heatmap + ggplot_abs_labs
}
}
ans <- list(cor_heatmap_list, cor_abs_heatmap_list)
return (ans)
}
load("SAVER_seurat_corr.RData")
# save(cormat_list, heatmap_list, saver_mat, sub_celltype, store_date, file = "SAVER_seurat_corr.RData")
# cormat_list <- make_cormat(saver_mat)
# heatmap_list <- make_cor_heatmap(saver_mat, cormat_list)
cor_pearson_mat <- cormat_list[[1]]; cor_spearman_mat <- cormat_list[[2]];
cor_kendall_mat <- cormat_list[[3]]; cor_hoeffd_mat <- cormat_list[[4]];
cor_blomqvist_mat <- cormat_list[[5]]; cor_MI_mat <- cormat_list[[6]];
cor_dist_mat <- cormat_list[[7]]; cor_XI_mat <- cormat_list[[8]];
1. Pearson’s correlation coefficient
- Pearson’s correlation is to measure linear dependency of data, X and Y
- \(-1 \leq \rho_{Pearson}(X, Y) \leq 1\)
- \(\rho_{Pearson}(X, Y) = \frac{\sum(x_i-\bar{x})(y_i -\bar{y})}{\sum(x_i-\bar{x})^2(y_i -\bar{y})^2}\)
cor_pearson_mat[1:5,1:5]
## TMSB4X MT-ND2 MT-ATP6 MT-CYB RPS25
## TMSB4X NA NA NA NA NA
## MT-ND2 -0.02940452 NA NA NA NA
## MT-ATP6 0.05764577 0.6005007 NA NA NA
## MT-CYB 0.01728402 0.6175195 0.6597476 NA NA
## RPS25 -0.45423356 -0.1214727 -0.1986036 -0.1419443 NA
quantile(cor_pearson_mat, na.rm = T)
## 0% 25% 50% 75% 100%
## -0.72766171 -0.08780746 0.02768556 0.14230112 0.96059754
quantile(abs(cor_pearson_mat), na.rm = T)
## 0% 25% 50% 75% 100%
## 3.489797e-07 5.478188e-02 1.171753e-01 2.024889e-01 9.605975e-01
# plot the smallest correlations
cor_pearson_vec <- sort(abs(cor_pearson_mat), decreasing = T)
plot(cor_pearson_vec)

#plot the high correlations
par(mfrow = c(2,2))
for(i in 1:4){
idx <- which(abs(cor_pearson_mat) == cor_pearson_vec[i], arr.ind = T)
idx1 <- idx[1]; idx2 <- idx[2]
plot(saver_mat[,idx1], saver_mat[,idx2], col = sub_cluster_labels, asp = T,
pch = 16, xlab = paste0(colnames(saver_mat)[idx1], ", (", idx1, ")"),
ylab = paste0(colnames(saver_mat)[idx2], ", (", idx2, ")"),
main = paste0("Correlation of ", round(cor_pearson_mat[idx1, idx2], 3)))
}

#plot the lowest correlations
par(mfrow = c(2,2))
for(i in 1:4){
idx <- which(abs(cor_pearson_mat) == rev(cor_pearson_vec)[i], arr.ind = T)
idx1 <- idx[1]; idx2 <- idx[2]
plot(saver_mat[,idx1], saver_mat[,idx2], col = sub_cluster_labels, asp = T,
pch = 16, xlab = paste0(colnames(saver_mat)[idx1], ", (", idx1, ")"),
ylab = paste0(colnames(saver_mat)[idx2], ", (", idx2, ")"),
main = paste0("Correlation of ", round(cor_pearson_mat[idx1, idx2], 3)))
}

Heatmap
heatmap_list[[1]][[1]]

2. Spearman’s correlation coefficient
- It captures the monotonic relationship between data, X and Y
- \(-1 \leq \rho_{Spearman}(X,Y) \leq 1\)
- \(\rho_{Spearman} = 1 - \frac{6\sum{d_i^2}}{n(n^2-1)}\) where \(d_i\) is the difference between the ranks of \(x_i\) and \(y_i\)
cor_spearman_mat[1:5,1:5]
## TMSB4X MT-ND2 MT-ATP6 MT-CYB RPS25
## TMSB4X NA NA NA NA NA
## MT-ND2 0.03179945 NA NA NA NA
## MT-ATP6 -0.05865502 0.5912240 NA NA NA
## MT-CYB -0.03411814 0.5749782 0.6511620 NA NA
## RPS25 -0.11125302 -0.1000161 -0.1893696 -0.109827 NA
quantile(cor_spearman_mat, na.rm = T)
## 0% 25% 50% 75% 100%
## -0.84618253 -0.07354277 0.02984100 0.13308323 0.81889279
quantile(abs(cor_spearman_mat), na.rm = T)
## 0% 25% 50% 75% 100%
## 1.613408e-06 4.973064e-02 1.053442e-01 1.810083e-01 8.461825e-01
# plot the smallest correlations
cor_spearman_vec <- sort(abs(cor_spearman_mat), decreasing = T)
plot(cor_spearman_vec)

#plot the high correlations
par(mfrow = c(2,2))
for(i in 1:4){
idx <- which(abs(cor_spearman_mat) == cor_spearman_vec[i], arr.ind = T)
idx1 <- idx[1]; idx2 <- idx[2]
plot(saver_mat[,idx1], saver_mat[,idx2], col = sub_cluster_labels, asp = T,
pch = 16, xlab = paste0(colnames(saver_mat)[idx1], ", (", idx1, ")"),
ylab = paste0(colnames(saver_mat)[idx2], ", (", idx2, ")"),
main = paste0("Correlation of ", round(cor_spearman_mat[idx1, idx2], 3)))
}

#plot the lowest correlations
par(mfrow = c(2,2))
for(i in 1:4){
idx <- which(abs(cor_spearman_mat) == rev(cor_spearman_vec)[i], arr.ind = T)
idx1 <- idx[1]; idx2 <- idx[2]
plot(saver_mat[,idx1], saver_mat[,idx2], col = sub_cluster_labels, asp = T,
pch = 16, xlab = paste0(colnames(saver_mat)[idx1], ", (", idx1, ")"),
ylab = paste0(colnames(saver_mat)[idx2], ", (", idx2, ")"),
main = paste0("Correlation of ", round(cor_spearman_mat[idx1, idx2], 3)))
}

Heatmap
heatmap_list[[1]][[2]]

3. Kendall’s Tau
- It is an alternative method to Spearman’s correlations, identifying monotonic relationships.
- \(-1 \leq \rho_{Kendall}(X,Y) \leq 1\)
- \(\rho_{Kendall}(X,Y) = \frac{\#\;concordant\;pairs - \#\;discordant \;pairs}{0.5n(n-1)}\)
cor_kendall_mat[1:5,1:5]
## TMSB4X MT-ND2 MT-ATP6 MT-CYB RPS25
## TMSB4X NA NA NA NA NA
## MT-ND2 0.01799150 NA NA NA NA
## MT-ATP6 -0.04475093 0.43496159 NA NA NA
## MT-CYB -0.02276791 0.40848512 0.4764382 NA NA
## RPS25 -0.07237012 -0.06244924 -0.1186702 -0.07486883 NA
quantile(cor_kendall_mat, na.rm = T)
## 0% 25% 50% 75% 100%
## -0.66004350 -0.05001927 0.02019916 0.08993943 0.63332765
quantile(abs(cor_kendall_mat), na.rm = T)
## 0% 25% 50% 75% 100%
## 3.962028e-06 3.370566e-02 7.132955e-02 1.228110e-01 6.600435e-01
# plot the smallest correlations
cor_kendall_vec <- sort(abs(cor_kendall_mat), decreasing = T)
plot(cor_kendall_vec)

#plot the high correlations
par(mfrow = c(2,2))
for(i in 1:4){
idx <- which(abs(cor_kendall_mat) == cor_kendall_vec[i], arr.ind = T)
idx1 <- idx[1]; idx2 <- idx[2]
plot(saver_mat[,idx1], saver_mat[,idx2], col = sub_cluster_labels, asp = T,
pch = 16, xlab = paste0(colnames(saver_mat)[idx1], ", (", idx1, ")"),
ylab = paste0(colnames(saver_mat)[idx2], ", (", idx2, ")"),
main = paste0("Correlation of ", round(cor_kendall_mat[idx1, idx2], 3)))
}

#plot the lowest correlations
par(mfrow = c(2,2))
for(i in 1:4){
idx <- which(abs(cor_kendall_mat) == rev(cor_kendall_vec)[i], arr.ind = T)
idx1 <- idx[1]; idx2 <- idx[2]
plot(saver_mat[,idx1], saver_mat[,idx2], col = sub_cluster_labels, asp = T,
pch = 16, xlab = paste0(colnames(saver_mat)[idx1], ", (", idx1, ")"),
ylab = paste0(colnames(saver_mat)[idx2], ", (", idx2, ")"),
main = paste0("Correlation of ", round(cor_kendall_mat[idx1, idx2], 3)))
}

Heatmap
heatmap_list[[1]][[3]]

4. Hoeffding’s D
- It tests the independence of data by calculating the distance between the product of the marginal distributions under the null hypothesis and the empirical bi-variate distribution.
- \(-1 \leq D(X,Y) \leq 1\)
- \(D(X,Y) = \frac{(n-2)(n-3)D_1+D_2-2(n-2)D_3}{n(n-1)(n-2)(n-3)(n-4)}\)
- \(D_1 = \sum_{i=1}^{n} Q_i(Q_i-1)\)
- \(D_2 = \sum_{i=1}^{n} (R_i-1)(R_i-2)(S_j-1)(S_j-2)\)
- \(D_3 = \sum_{i=1}^{n} (R_i-2)(S_i-2)Q_i\)
cor_hoeffd_mat[1:5,1:5]
## TMSB4X MT-ND2 MT-ATP6 MT-CYB RPS25
## TMSB4X NA NA NA NA NA
## MT-ND2 0.009290504 NA NA NA NA
## MT-ATP6 0.011479529 0.135400141 NA NA NA
## MT-CYB 0.006091448 0.109296103 0.15678142 NA NA
## RPS25 0.008666603 0.007664955 0.01265554 0.005907249 NA
quantile(cor_hoeffd_mat, na.rm = T)
## 0% 25% 50% 75% 100%
## 0.0006317592 0.0052005399 0.0082638904 0.0145437026 0.3361010002
# plot the smallest correlations
cor_hoeffd_vec <- sort(abs(cor_hoeffd_mat), decreasing = T)
plot(cor_hoeffd_vec)

#plot the high correlations
par(mfrow = c(2,2))
for(i in 1:4){
idx <- which(abs(cor_hoeffd_mat) == cor_hoeffd_vec[i], arr.ind = T)
idx1 <- idx[1]; idx2 <- idx[2]
plot(saver_mat[,idx1], saver_mat[,idx2], col = sub_cluster_labels, asp = T,
pch = 16, xlab = paste0(colnames(saver_mat)[idx1], ", (", idx1, ")"),
ylab = paste0(colnames(saver_mat)[idx2], ", (", idx2, ")"),
main = paste0("Correlation of ", round(cor_hoeffd_mat[idx1, idx2], 3)))
}

#plot the lowest correlations
par(mfrow = c(2,2))
for(i in 1:4){
idx <- which(abs(cor_hoeffd_mat) == rev(cor_hoeffd_vec)[i], arr.ind = T)
idx1 <- idx[1]; idx2 <- idx[2]
plot(saver_mat[,idx1], saver_mat[,idx2], col = sub_cluster_labels, asp = T,
pch = 16, xlab = paste0(colnames(saver_mat)[idx1], ", (", idx1, ")"),
ylab = paste0(colnames(saver_mat)[idx2], ", (", idx2, ")"),
main = paste0("Correlation of ", round(cor_hoeffd_mat[idx1, idx2], 3)))
}

Heatmap
heatmap_list[[1]][[4]]

5. Blomqvist’s Beta
- It measures dependency between variables by constructing a two-way contingency table with the medians of each margin as cutting points.
- \(0 \leq \beta \leq 1\)
- \(\beta_n = \frac{n_1-n_2}{n_1+n_2} = \frac{2n_1}{n_1+n_2} - 1\)
- \(\beta = P\{(X-\tilde{x})(Y-\tilde{y})>0\} - P\{(X-\tilde{x})(Y-\tilde{y}) < 0\}\)
cor_blomqvist_mat[1:5,1:5]
## [,1] [,2] [,3] [,4] [,5]
## [1,] NA NA NA NA NA
## [2,] 1 NA NA NA NA
## [3,] 1 1 NA NA NA
## [4,] 1 1 1 NA NA
## [5,] 1 1 1 1 NA
quantile(cor_blomqvist_mat, na.rm = T)
## 0% 25% 50% 75% 100%
## 0.9881188 1.0000000 1.0000000 1.0000000 1.0000000
quantile(abs(cor_blomqvist_mat), na.rm = T)
## 0% 25% 50% 75% 100%
## 0.9881188 1.0000000 1.0000000 1.0000000 1.0000000
# plot the smallest correlations
cor_blomqvist_vec <- sort(abs(cor_blomqvist_mat), decreasing = T)
plot(cor_blomqvist_vec)

#plot the high correlations
par(mfrow = c(2,2))
for(i in 1:4){
idx <- which(unique(abs(cor_blomqvist_mat)) == unique(cor_blomqvist_vec)[i], arr.ind = T)
idx1 <- idx[1,1]; idx2 <- idx[1,2]
plot(saver_mat[,idx1], saver_mat[,idx2], col = sub_cluster_labels, asp = T,
pch = 16, xlab = paste0(colnames(saver_mat)[idx1], ", (", idx1, ")"),
ylab = paste0(colnames(saver_mat)[idx2], ", (", idx2, ")"),
main = paste0("Correlation of ", round(cor_blomqvist_mat[idx1, idx2], 3)))
}

#plot the lowest correlations
par(mfrow = c(2,2))
for(i in 1:4){
idx <- which(abs(cor_blomqvist_mat) == rev(cor_blomqvist_vec)[i], arr.ind = T)
idx1 <- idx[1,1]; idx2 <- idx[1,2]
plot(saver_mat[,idx1], saver_mat[,idx2], col = sub_cluster_labels, asp = T,
pch = 16, xlab = paste0(colnames(saver_mat)[idx1], ", (", idx1, ")"),
ylab = paste0(colnames(saver_mat)[idx2], ", (", idx2, ")"),
main = paste0("Correlation of ", round(cor_blomqvist_mat[idx1, idx2], 3)))
}

Heatmap
heatmap_list[[1]][[5]]
## NULL
6. Normalized Mutual Information
- It measures how much one random variable gives information about the other. For example, High mutual information indicates a large reduction in uncertainty.
- \(0 \leq MI(X,Y) \leq 1\), as it is normalized.
- $MI(X,Y) = f_{X,Y} (x,y) log_2 ; dxdy $
- \(MI(X,Y) = \sum \sum p_{X,Y} (x,y) log \frac{p_{X,Y} (x,y)}{P_X(x)P_Y(y)}\)
cor_MI_mat[1:5,1:5]
## [,1] [,2] [,3] [,4] [,5]
## [1,] NA NA NA NA NA
## [2,] 0.5900765 NA NA NA NA
## [3,] 0.6587018 0.2397348 NA NA NA
## [4,] 0.7150763 0.2277529 0.3148607 NA NA
## [5,] 0.9971672 0.1872224 0.0834614 0.1378969 NA
quantile(cor_MI_mat, na.rm = T)
## 0% 25% 50% 75% 100%
## 0.02008612 0.11215506 0.13350667 0.15808349 1.00000000
# plot the smallest correlations
cor_MI_vec <- sort(abs(cor_MI_mat), decreasing = T)
plot(cor_MI_vec)

#plot the high correlations
par(mfrow = c(2,2))
for(i in 1:4){
idx <- which(unique(abs(cor_MI_mat)) == unique(cor_MI_vec)[i], arr.ind = T)
idx1 <- idx[1,1]; idx2 <- idx[1,2]
plot(saver_mat[,idx1], saver_mat[,idx2], col = sub_cluster_labels, asp = T,
pch = 16, xlab = paste0(colnames(saver_mat)[idx1], ", (", idx1, ")"),
ylab = paste0(colnames(saver_mat)[idx2], ", (", idx2, ")"),
main = paste0("Correlation of ", round(cor_MI_mat[idx1, idx2], 3)))
}

#plot the lowest correlations
par(mfrow = c(2,2))
for(i in 1:4){
idx <- which(abs(cor_MI_mat) == rev(cor_MI_vec)[i], arr.ind = T)
idx1 <- idx[1]; idx2 <- idx[2]
plot(saver_mat[,idx1], saver_mat[,idx2], col = sub_cluster_labels, asp = T,
pch = 16, xlab = paste0(colnames(saver_mat)[idx1], ", (", idx1, ")"),
ylab = paste0(colnames(saver_mat)[idx2], ", (", idx2, ")"),
main = paste0("Correlation of ", round(cor_MI_mat[idx1, idx2], 3)))
}

Heatmap
heatmap_list[[1]][[6]]
